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Abstract 

We study the dynamics of a SIR epidemic model with nonlinear incidence 
rate, vertical transmission vaccination for the newborns and the capac¬ 
ity of treatment, that takes into account the limitedness of the medical 
resources and the efficiency of the supply of available medical resources. 

Under some conditions we prove the existence of backward bifurcation, 
the stability and the direction of Hopf bifurcation. We also explore how 
the mechanism of backward bifurcation affects the control of the infectious 
disease. Numerical simulations are presented to illustrate the theoretical 
findings. 

Keywords: Local Stability, Hopf Bifurcation, Global Stability, Backward Bi¬ 
furcation 


1 Introduction 

Mathematical models that describe the dynamics of infectious diseases in com¬ 
munities, regions and countries can contribute to have better approaches in the 
disease control in epidemiology. Researchers always look for thresholds, equilib¬ 
ria, periodic solutions, persistence and eradication of the disease. For classical 
disease transmission models, it is common to have one endemic equilibrium and 
that the basic reproduction number tells us that a disease is persistent if is 
greater than Rand dies out if is less than 1. This kind of behaviour associates 
to forward bifurcation. However, there are epidemic models with multiple en¬ 
demic equilibria [1, 2, 3, 4], within these models it can happen that a stable 
endemic equilibrium coexist with a disease free equilibrium, this phenomenon 
is called backward bifurcation [5]. 

In order to prevent and control the spread of infectious diseases like, measles, 
tuberculosis and influenza, treatment is an important and effective method. In 
classical epidemic models, the treatment rate of the infectious is assumed to be 
proportional to the number of the infective individuals [6]. Therefore we need to 
investigate how the application of treatment affects the dynamical behaviour of 
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these diseases. In that direction in [7], Wang and Ruan, considered the removal 
rate 


T(I) = 


k, if I > 0 
0 , if I = 0 . 


In the following model 


— = A-dS- XSI, 
dt 

±L=\SI-(d + 'y)I-T(I), 

f = 7 / + r ( /)-dfl. 


where S, I, and R denote the numbers of the susceptible, infective and recovered 
individuals at time t , respectively. The authors study the stability of equilibria 
and prove the model exhibits Bogdanov-Takens bifurcation, Hopf bifurcation 
and Homoclinic bifurcation. In [15], the authors introduce a saturated treatment 


T(I) 


PI 

1 H- olI 


A related work is [13], [16]. 

Hu, Ma and Ruan [ 8 ] studied the model 


dS 

dt 

dJ 

dt 

dR 

dt 


=bm(S + R) — 


/3SI 
l + al 


-bS+pSI 


psi 


+ (qS - 6 - 7 )/ - T(I) 


1 -b OlI 

= 7 / - bR + bm'(S + R) + T(I) 


a) 


the basic assumptions for the model ( 1 ) are, the total population size at time 
t is denoted by N = S + I + R. The newborns of S and R are susceptible 
individuals, and the newborns of I who are not vertically infected are also 
susceptible individuals, b denotes the death rate and birth rate of susceptible 
and recovered individuals, 6 denotes the death rate and birth rate of infective 
individuals, 7 is the natural recovery rate of infective individuals, q (q < 1 ) is 
the vertical transmission rate, and note p = 1 — q, then 0 < p < 1. Fraction 
m! of all newborns with mothers in the susceptible and recovered classes are 
vaccinated and appeared in the recovered class, while the remaining fraction, 
m = 1 — m', appears in the susceptible class, the incidence rate is described by 
a nonlinear function /3S7/(1 + al), where p is a positive constant describing 
the infection rate and a is a nonnegative constant. The treatment rate of the 
disease is 


T(I) = 


hi, if 0 < / < / 0 , 
u = klo, if I > Iq 


where Iq is the infective level at which the healthcare systems reaches capacity. 
In this work we will extend model (1) introducing the treatment rate - 21 


1 +a 2 I ’ 
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where a 2 , ft 2 > 0 , obtaining the following model 


dR 

dt 

dJ 

dt 

dR 

dt 


=bm(S + R) — 
/3SI 


PSI 

l + al 


l + aR ^ 5 - 5 -' 1 ' 1 
= 7 1 — bR + bm'(S + R) - 


bS + P 5I 

_ fal 
1 + a 2 I 
021 

1 + a 2 f 


( 2 ) 


Because ^ = 0, the total number of population N is constant. For convenience, 
it is assumed that N = S + I + R=1. By using S + R = 1 — I, the first 
two equations of (2) do not contain the variable R. Therefore, system (2) is 
equivalent to the following 2 -dimensional system: 


dS 

dt 

dJ 

dt 


psi 

1 + al 


bS + bm( 1 — I) + pSI 


psi 

1 H- cxI 


pSI — 7 / 


P 2 I 

1 + a 2 1 


( 3 ) 


The parameters in the model are described below: 


• S, /, R are the normalized susceptible, infected, and recovered population, 
respectively, therefore it follows that S,I,R< 1. 

• b is a positive number representing the birth and death rate of susceptible 
and recovered population. 

• S is a positive number representing the birth and death rate of infected 
population. 

• 7 is a positive number giving the natural recovery rate of infected popu¬ 
lation. 

• q is positive (q < 1) representing the vertical transmission rate (disease 
transmission from mother to son before or during birth). It is assumed 
that descendents of the susceptible and recovered classes belong to the 
susceptible class, in the same way to the fraction of the newborns of the 
infected class not affected by vertical transmission. 

• p = 1 — q therefore 0 < p < 1. 

• m! is positive and it is the fraction of vaccinated newborns from susceptible 
and recovered mothers and therefore belong to the recovered class, m = 
1 — ml > 0 is the rest of newborns, which belong to the susceptible class. 

• P is positive, representing the infection rate, cc is a positive saturation 
constant (In the model the incidence rate is given by the nonlinear function 

MI \ 

1 -\-OtI ' 
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• the treatment function, satisfying lim/-^ where 

Ot2, P2 > 0 . 

We note that if a 2 = 0 the treatment becomes bilinear, case considered in 
[ 8 ], whereas if /?2 = 0 treatment is null, not being of interest here. Therefore we 
will assume /3 2 ,a 2 > 0 . 

The paper in distributed as follows: in section 2 we compute the equilibria 
points and determine the conditions of its existence (as real values) and posi¬ 
tivity, in section 3 we analyze the stability of the disease free equilibrium and 
endemic equilibria points in terms of value of TZq and the parameters of treat¬ 
ment function. Section 4 is dedicated to study Hopf bifurcation of the endemic 
equilibria points and section 5 shows discussion of all our results and we give 
some control measures that could be effective to eradicate the disease in each 
case. 

Following [ 8 ] we define 


IZo := 


/3m 

(32 + pS + 7 


When /3 2 = 0, TZq reduces to 


n 


* _ 

0 — 


/3 m 

pS + 7 ’ 


(4) 

(5) 


which is the basic reproduction number of model (3) without treatment. 

Lemma 1 . Given the initial conditions S( 0) = So > 0, 1(0) = Iq > 0, then the 
solution of (3) satisfies 5(f),/(f) > 0 Vt > 0 and S(t) + I (t) < 1. 

Proof. Take the solution 5(f),/(f) satisfying the initial conditions 5(0) = So > 
0,/(0) = Jo > 0. Assume that the solution is not always positive, i.e., there 
exists a to such that 5(to) < 0 or /(to) < 0. By Bolzano’s theorem there 
exists a ti € (0, to] such that 5(ii) = 0 or /(ti) = 0 , which can be written as 
5(ti)/(<i) = 0 for some ti € (0, to]. Let 


t 2 = min {ti,S(ti)I(ti) = 0 }. ( 6 ) 

Assume first that 5(t 2 ) = 0, then dS ^ > 0 implying that 5 is increasing at t = 
t 2 . Hence 5(f) is negative for values of f < f 2 near f 2 , a contradiction. Therefore 
5(f) > 0 Vt > 0 and we must have /(t 2 ) = 0, implying = 0. Note that if 
for some f > 0 J(t) = 0, then = 0. Then any solution with 1(0) — Iq = 0 
will satisfy /(f) = 0 Vt > 0. By uniqueness of solutions this fact implies that if 
1(0) = Iq > 0, then /(f) will remain positive for all t > 0. Therefore /(f 2 ) = 0 
leads to a contradiction. Hence both 5 and I are nonnegative for all t > 0. 
Finally, adding both derivatives of 5(f) and I(t) we get: 


d(S + I) 
dt. 


—bS + bm — bml — 7 / 


P2I 

1 + a 2 I 


(7) 
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Being S, I > 0, if S+I — 1 then 0 <iS< 1,0</<1. Analyzing the expression 
—bS + bm — bml, 

—bS + bm- — bml = b{m — ml — S) = b{m — ml — 1 + 1) = b(m — 1 + 1(1 — to)). 

Note that by the definition of the model parameters, 1 — m = m' > 0. Knowing 
that / < 1, then 


/(I — m) < 1 — to => 1(1 — to) + to — 1 < 0 . ( 8 ) 

Therefore — bS + bm. — bml < 0. Hence < 0 and S + I is non increasing 

along the line S +1 = 1, implying that S +1 < 1. Note also that S + I cannot 
be grater than 1, otherwise from R = 1 — (S + /), R would be negative, a 
nonsense. □ 


2 Existence and positivity of equilibria 

Assume that system (3) has a constant solution (Sq,Iq), then: 


PSqIq 
1 + al 0 
PS 0 I 0 


— bS 0 + bm( 1 — I 0 ) + pSI 0 = 0 




From (9) we obtain 


Sn = 


(1 + al 0 )(bm(l - ip) + pSI 0 ) 

(3Io + b(l + alo) 


And we get from (10): 

(3S 0 


1 + CxIq 
=> I 0 = 0 or 


— p8 — 7 — 


& 


1 + « 2^0 


= 0 


PSo @2 n 


(9) 

( 10 ) 

( 11 ) 


( 12 ) 


If Iq = 0 then Sq = to, obtaining in that way the disease-free equilibrium 
E = (?n, 0). 

Theorem 2. System (3) has a positive disease-free equilibrium E = (to, 0). 

In order to obtain positive solutions of system 3 if Iq ^ 0 then: 


— p5 — 7 — 


02 


ps 0 

1 + alo 

=> So = 1 + ^ Ia (pS + 7 


1 + CX2I0 


= 0 


02 


1 + 02/0 


(13) 
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We obtain the following quadratic equation: 


AIq + BIq + C = 0. 

(14) 

Or 


Iq + ( B/A)Iq + C/A = 0, 

(15) 


where the coefficients are given by: 


A = a 2 (P( 7 + bm ) + ab(pS + 7)) > 0 

B = /?( 7 + p 2 + bm(l — a 2 )) + ba(jp5 + 7 + / 3 2 ) + ba 2 (p6 + 7) 

= /?(7 + /?2 + bm — bma 2 ) + ba( 1 — TZ 0 )(pS + 7 + /? 2 ) + /3rn&a + ba 2 (pS + 7) 
C = b(pS + 7 + /3 2 - pm) = b(pS + 7 + /3 2 )(1 - U 0 ). (16) 


Its roots are: 


h 

h 


B - sjB- - 4AC 
2 A 

-B + VB 2 - 4AC 
2A 


And using these values in (13) we obtain its respective values 


(17) 


Si 


S 2 


1 + ah 

1 + al 2 


P 


( pi+7+ iT 3 j kr) 


(18) 


Then our candidate for endemic equilibria are Ei = (Si,h), E 2 = ( S 2l I 2 ). 

Note that C = 0 if and only if 77 q = 1, C > 0 if and only if 77o < 1, and 
C < 0 if and only if 77o > 1 . 

For TZq > 1 we define the following sets: 


A\ = {(P 2 , a 2 ) : p 2 > 0,0 < a 2 < a®, 

A 2 = {(P 2 ,a 2 ) : p 2 > g(a 2 ),a 2 > > 0}, 

A 3 = {(P 2 , a 2 ) : 0 < p 2 < g{a 2 ),a 2 > > 0 }. (19) 


Where 


g(a 2) 


q — P(mba + 7 + bm) 

b(p 5 + 7 — Pm) 


1 

P 


( ba 2 (pS + 7 — pm) + P( 7 + bm + mba)). 


Define : 
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Figure 1: Location of the sets A ll A 2 ,A 3 in the plane a 2 — /3 2 , for 7 = 0.01, /3 = 
0.2, b = 0.2, m = 0.3, p = 0.02,5 = 0.1 


Pi = 1 


1 


P,T = 1 - 


ba(p6 + 7 + /3 2 ) 
1 


[/? (7 + P 2 + bin — bma 2 ) + (3mba + ba 2 (p6 + 7 )] 


ba 2 (p6 + 7 + /3 2 ) 

\J—/3a(bma + /3 2 + 7 + bm — a 2 bm ) + /3a 2 ( r y + bm) — \J a 2 (/J 7 + (3bm + abp5 + ab 7 ) 

( 20 ) 


Figure 1 shows the location of these sets. 

Theorem 3. If TZq > 1 the system (3) has a unique (positive) endemic equilib¬ 
rium E 2 . 

Proof. If TZq > 1 then C < 0, then using Routh Hurwitz criterion for n = 2, the 
quadratic equation has two real roots with different sign, /1 and I 2l where I\ < 
I 2 . Hence there exists a unique positive endemic equilibrium E 2 = (S 2 ,I 2 ). □ 


Theorem 4. Let 0 TZq 1 . For system (3), if TZq ^ 1 then there are no 
positive endemic equilibria. Otherwise, if TZq > 1 the following propositions 
hold: 


1. If TZq = 1 and (/3 2 ,a 2 ) € A 3 the system (3) has a unique positive endemic 
equilibrium E 2 = ( S 2l I 2 ), where 

,^-B/A. & = l^( trf + 7+n ^). 
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2. If max{Pi, Rq } < 72q < 1 and ( 02 , (* 2 ) G A3, the system (3) has a pair of 
positive endemic equilibria Pi, P 2 . 

3. If 1 > IZo = Rq > Pi and ( 02 , 02) G A 3 , the system (3) has a unique 
positive endemic equilibrium E\ = P 2 . 

4■ If 1 > 72.0 = Pi an d ( 02 , <* 2 ) G A 3 , f/ie system (3) has no positive endemic 
equilibria. 

5. If 0 < 72o < 1 and ( 02 , 012 ) G Ai U ^ 2 , f/ie system (3) has no positive 
endemic equilibria. 

6 . If ( 02 , 02 ) G A 3 and 0 < 72o < max(7?g",Pi) < 1, then there are no 
positive endemic equilibria. 


Proof. If 0 < 72-0 < 1, then C > 0, so the roots of the equation AI 2 + BI+C = 0 
are not real with different sign, but real with equal signs, complex conjugate 
or some of them are zero. If endemic equilibria exist and are positive, it is 
necessary that B < 0. After some calculations we can see that: 


P<0<S=> 72 0 >1 + 


0 (l + 02 + bm — bma 2 ) + 0 mba + ba 2 (pS + 7 ) 
ba(p 6 + 7 + 0 2 ) 



From the assumption that 72o < 1 then Pi < 1, hence the expression 0("f + 02 + 
bm — bm. 012 ) + 0 mba + ba 2 (pS + 7 ) must be negative, this happens if and only if 


02 < 


1 

0 


(ba 2 (p 5 + 7 — 0 m) + 0 (j + bm + mba)) = < 7 ( 0 : 2 ). 


( 22 ) 


If 72 q < 1 then —^( 6 o 2 (p^ + 7 — 0m) + 0(j + bm + mba)) < 0 and it is not 
possible to find a value of 02 fulfilling the previous inequality, therefore there 
are no positive endemic equilibria. 

Now, if 72 q > 1 we have that. 


1. If 72o = 1 then C = 0 and the equation (14) is transformed into 

AIq + Bio = 0, 


(23) 


with A > 0. Its roots are I\ = 0 and I 2 = —B/A, and there exists a 
unique endemic equilibrium that is positive if and only if B < 0, that is 
given by E 2 = (<S 2 ,/ 2 ), where 


I2 

S2 


-B/A 
1 + 0/2 
~ 1 ~ 


pS + 7 + 


02 \ 
1 + a 2 / 2 ) 


(24) 


Note that if o 2 > 0 ° and 72g > 1 then g(a 2 ) > 0. 

Hence A 3 is nonempty and its elements satisfy B < 0, therefore if ( 02 , 02 ) G 
A 3 there exists a unique positive endemic equilibrium p 2 - 



2. If TZ 0 < 1 then C > 0 and the roots of the quadratic equation for I 0 must 
be real of equal sign or complex conjugate. By the previous part we know 
that if ( 02 , 012 ) G A 3 then Pi < 1, moreover if 1Z 0 > Pi then B < 0 and 
therefore both roots must have positive real part. Finally, to assure that 
equilibria are both real, we demand that A > 0 . Computing A: 

A = B 2 - 4AC 

= A 2 K 2 0 + B 2 K 0 +C 2 = A(7Z 0 ), (25) 

where: 

A 2 = a 2 b 2 (p 6 + 7 + 0 2 ) 2 (26) 

B 2 = —2 \0 (7 + 02 + bm (1 — as)) + a b (pS + 7 + 0s) + 0 mba 
+ ba 2 (pS + y)]ab (pS + 7 + 0s) + 4a 2 (0 (7 + bm) a b (pS + 7)) 
b(pS + y + 0 2 ) (27) 

C 2 = (0 (7 + 0s + bm (1 - as)) +ab(p5 + 7 + 0 2 ) + 0mba + ba 2 (pS + 7 )) 2 
— 4q 2 (0 (7 + bm) + ab (pS + 7)) b (pS + 7 + 0s) ■ (28) 


The previous expression is a quadratic function of TZq. To establish the 
region where A > 0, it is necessary to know how the roots of A(77. 0 ) 
behave. The discriminant of the quadratic function A(77o) is 

A 2 = —16 a 2 b 2 0 (pS + 7 + 0s) 2 (01 + 0 bm + a bp5 + a by) 

(a(abm + 0 2 + 7 + bm) — a 2 (y + bm- + abm)). (29) 


If we assume that A 2 < 0, then a 2 < and in this case we 

have that: 


7 + 02 + bm — bma 2 + bma > 


20 2 abm + (7 + 6 m) (7 + 02 + bm + bma) 


7 + bm + abm 


> 0 . 


(30) 


So we get that P\ > 1 > 77o, which is a contradiction with the assumption 
in this part, therefore A 2 > 0 and in consequence A(77o) has two real 
roots, 


R 0 — 


—B 2 — \J A 2 


2A 2 


= 1 - 


ba' 2 (pS + 7 + 0 2 ) 

+ \J a 2 ( 0 (y + bm) + ab(pS + 7 ))] 2 , 


[\/— 0 (a(bma + 0 2 + 7 + bm — bma 2 ) — a 2 (7 + bm)) 


+ _ —B 2 + sfKs 

rir\ — 


2Ao 


1 - 


ba 2 (p 6 + 7 + 0 2 ) 

- \Ja 2 (0(y + bm) + ab(pS + 7))] 2 


[\/— 0 (a(bma + 0 2 + 7 + bm- — bma 2 ) — a 2 (y + bm)) 

(31) 
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Note that due to the positivity of A 2 and (29), we have that 

—/ 3 (a(bma + /?2 + 7 + bm — bma 2 ) — a 2 (7 + bm)) 

is positive, allowing its roots to be well defined. Analyzing the derivative 
of A (Rq) we have that 

A'(R+) = > 0 and A'(Rq ) = < 0 , 

moreover Rq < R q making A positive for Rq > Rq or Rq < Rq . Never¬ 
theless 

Rq = 1 + , , . 1 —T 77 t (^(7 + P 2 + bm - bma 2 + bma)) - e, 
ba(pd + 7 + p 2 ) 


while 


Pi = l + 


1 

ba(pS + 7 + /? 2 ) 


(/3 (7 + /3 2 + bm — bma .2 + bma)) + e 2 , 


with e, e 2 > 0, making R g < Pi < R 0 . Therefore for 1Z 0 > max(Pi, Rq), 
we have that there exists two positive endemic equilibria E \ 1 E 2 , proving 
this part. 

3. If (p 2 ,ct 2 ) € A 3 then Pi < 1. If 1 > Pq > Pi, then we have that B < 0 
and C > 0, therefore we have a pair of roots of the quadratic for I with 
positive real part. In the previous part it was proven that for Pi < 1 the 
discriminant A 2 > 0 and both roots Rq , Rq are real and less than one. If 
1Z 0 = Rq then A = 0 and both roots are fused in one Ii = —B/2A = J 2 . 
Therefore we have a unique positive endemic equilibrium E 1 = E 2 - 

4. If (/3 2 ,a 2 ) € A 3 then Pi <1. If Pq = P| < 1 then C > 0, implying 
that the roots are complex conjugate or real of the same sign. Being 
TZ 0 = P 1 then B = 0, implying that both roots have real part equal to 
zero, therefore there are no positive endemic equilibria. 

5. If 0 < Rq < 1 and (/3 2 ,a 2 ) € Ai U A 2 then Pi > 1, therefore Rq < Pi 
, B > 0, and C > 0. Hence there are two roots with real part zero or 
negative, which are not positive equilibria. 

6 . If (P 2 , (* 2 ) G A 3 we have that Pi < 1 and the roots of the discriminant 
Rq , Rq are real, in addition that Rq < Pi and Rq < 1 by definition of 
this case. If 0 < Rq < max{PQ ,Pi} < 1, then C > 0 and the roots / 2 ,/3 
are complex conjugate or real with the same sign. If Rq < Pi then B > 0, 
and the roots have negative real part, so there are not positive endemic 
equilibria. If 0 < Rq < Rq and Rq > Rq , then A < 0 and the roots 
are complex conjugate, therefore there is not real endemic equilibria. If 
0 < Rq < Rq and Rq < Rq < Pi, then it reduces to the first case in 
which there are not positive endemic equilibria. 
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Figure 2: Graph of TZq versus the values I of equilibria. Parameter values used 
are a = 0.4, 0.2 = 3.8,/3 = 0.2,6 = 0.2,7 = 0.01,6 = 0.01, p = 0.02, to = 0.1. 
In this example /?2 varies from 0 to 0.025, therefore TZq varies between 0.5682 
and 1.9682. g(a.2 ) = —0.0017 and a® = 3.8776, so (^ 2 , 0 : 2 ) £ A-| U A 2 . Forward 
bifurcation can be observed in TZq = 1. 


Theorem 4 gives us a complete scenario of the existence of endemic equilibria. 
When TZq < 1 we have that TZq < 1, it follows from the fact that TZ 0 < TZq 
whenever @2 > 0; then system 3 has only a disease free equilibrium and no 
endemic equilibria. 

Otherwise, when TZq > 1 . If (/? 2 ,« 2 ) £ A x U A 2 then we have no endemic 
equilibria for 0 < TZq < 1 and a unique endemic equilibria E 2 when TZq > 1 , so 
there exists a forward bifurcation in TZq = 1 from the disease free equilibrium 
to E 2 (see figure 2 ). If (/3 2 , 0 : 2 ) £ A 3 there exist two positive endemic equilibria 
whenever max{Pi,.R^"} < 77-o < 1 ( Pi and TZq depend on /3 2 ), we can observe 
the backward bifurcation of the equilibrium E to two endemic equilibria (see 
figure 3 ). 

As an immediate consequence of the previous theorem we have that if TZq > 1 
there exists a unique positive endemic equilibrium, while if TZq < 1 and the 
conditions of the second part are fulfilled, there exist two positive endemic 
equilibria. Hence we have the following corollary: □ 

Corollary 5. IfTZ 0 = 1, TZq > 1 and (/?, 02 ) £ Aq, system (3) has a backward 
bifurcation of the disease-free equilibrium E. 

Proof. First we note that if (/3 2 5 0 - 2 ) £ ^3 then Rq is real less than one and Pi < 
1, therefore we can find a neighborhood of points in the interval (max{P^", Pi}, 1). 
By case 2 of theorem, if TZq lies in this neighborhood there exist two positive 
endemic equilibria P-|, E 2 : for TZq = 1 there exists a unique positive endemic 


11 


























Figure 3: Graph of TZq versus the values I of equilibria. In this example /3 2 varies 
from 0 to 0.025 therefore IZq varies between 0.5682 and 1.9682. Parameter values 
used are a = 0.4, a-i = 16, /3 = 0.2 ,6 = 0.2 ,7 = 0.01, (5 = 0.01, p = 0.02, m = 0.1. 
g(a 2 ) = 0.1188 and a 2 = 3.8776, so (fi 2 ,ct 2 ) £ A 3 . Backward bifurcation can 
be observed in TZq = 1 and the existence of two positive endemic equilibria 
whenever max{P!, Rq } < lZ a < 1 . 


equilibrium E 2 , while the other endemic equilibrium becomes zero. Finally for 
TZq > 1 there exists a unique positive endemic equilibrium as the zero ’’endemic” 
equilibrium becomes negative. O 


3 Characteristic Equation and Stability 


The characteristic equation of the linearization of system (3) in the equilibrium 


(So,/ 0 ) is given by: 

det(DF - 

XI), 

(32) 

where 

f dfi 

dfi\ 



DF=(§g 

\ OS 

&)■ 

dl / 

(33) 


Matrix is evaluated in the equilibrium (Sq,Iq). Functions /i ,/2 are the follow¬ 
ing: 


h = 


/3SI 

l + al 


bS + &7tt(1 — /) + pSI 


h 


pSI 

1 + al 


p5I — 7 / 


P2I 

1 + a 2 I 


(34) 

(35) 


Computing the matrix DF we obtain: 
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DF(S,I) 


(36) 


-01 



-0 s 


0S 


(1 + oT)' 


— bm. + pS 


02 


(1 + al) 2 p6 7 (1 + a 2 iy 


3.1 Stability of disease free equilibrium 

For the disease free equilibrium E = (to, 0) the Jacobian matrix is: 


DF(m, 0) = 


—b —0m — bm + pS 
0 0 in — p5 — 7 — 02 


Theorem 6 . IflZo < 1 then the equilibrium E = (to, 0) of model (3) is locally 
asymptotically stable, while iflZo > 1 then it is unstable. 

Proof. The characteristic equation for the equilibrium E is given by 


P(A) =det(DF(TO,0)-A/ 2x2 ) 

_ , f—b — A — 0m — bm + pS \ 

c ° y 0 0m — pS — 7 — 02 — Ay 

= (-b - \){0m -pS - 7 - 02 - A). (37) 

The equation (37) has two real roots Ai = —b and A 2 = 0m — pS — 7 — /3 2 . By 
Hartman-Grobman’s theorem, if the roots of (37) have non-zero real part then 
the solutions of system (3) and its linearization are qualitatively equivalent. If 
both roots have negative real part then the equilibrium E is locally asymptot¬ 
ically stable, whilst if any of the roots has positive real part the equilibrium is 
unstable. Clearly Ai < 0, but A 2 < 0 if and only if 


0m — pS - 7 < 02, 


if and only if TZq <1. □ 

According to the previous theorem and theorem 4 we obtain the following 
result for the global stability of equilibrium E : 

Theorem 7. If 0 < IZo < 1 and one of the following conditions holds: 

• n* 0 < 1. 

• 1Z 0 = Pi and (02, 0 : 2 ) S A 3 . 

• ( 02 , £*2) £ AiU A 2 . 

• (02,ct2) £ A 3 and 0 < TZq < max{i?Q , Pi}. 

Then equilibrium E of system (3) is globally asymptoticaly stable. 
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Figure 4: Global stability of equilibrium E. 


Proof. If 0 < IZq < 1 then by theorem 6 the equilibrium E is locally asymptoti¬ 
cally stable. If any of the given conditions holds then by theorem 4 there are no 
endemic equilibria in the region D = {S(t), I it) >0 Vf > 0, S(t) + I(t) < 1}, 
which it was proven to be positively invariant in theorem 1. By [12] (page 245) 
any solution of (3) starting in D must approach either an equilibrium or a closed 
orbit in D. By [14] (theorem 3.41) if the solution path approaches a closed or¬ 
bit, then this closed orbit must enclose an equilibrium. Nevertheless, the only 
equilibrium existing in D is E and it is located in the boundary of D , therefore 
there is no closed orbit enclosing it, totally contained in D. Hence any solution 
of system (3) with initial conditions in D must approach the point E as t tends 
to infinity. □ 

Example 8. Take the following values for the parameters: a = 0.4, a .2 = 
10, /3 = 0.2, b = 0.2, 7 = 0.01, 6 = 0.01, p = 0.02, m = 0.3, /?2 = 0.1. Equilibrium 
E = (0.3, 0) , IZq = 0.5445 < 1. By theorem 6, E is locally asymptotically stable, 
= 7-42 < «2 and g(a 2 ) = —0.1864 < /3 2 , therefore (/? 2 , < 22 ) £ A 2 . By theorem 
4 there are no positive endemic equilibria. Finally by theorem 7 we have that E 
is globally stable. See figure 4- 

Theorem 9. IfTZo = 1 and (3 2 7 ^ g(a 2 ) then equilibrium E is a saddle point. 
Moreover, if (/? 2 , ex 2 ) £ A\ U A 2 the region D is contained in the stable manifold 
of E. 

Proof. If IZo = 1 one of the eigenvalues of the Jacobian matrix of the system is 
zero, hence we cannot apply Hartman-Grobman’s theorem. In order to establish 
the stability of equilibrium E we apply central manifold theory. Making the 


14 





























change of variables, S = S — m, I = I , we obtain the equivalent system 


dS 

dt 

dl 

dt 


/3{S + m)I 

1 A Oil 

(3(S + m)I _ 
1 -t~ Cxi 


— bS — bml + p5i 

hi 


p5I — 7 / — 


1 + CX 2 I 


( 38 ) 


Because I — I we ignore the hat and use only I. This new system has an 
equilibrium in E = (0,0) and its Jacobian matrix in that point is 

DF(m, 0)=(- 6 -Pm-bm+pSy (39) 

Using change of variables S = u— A+ beta ^+ bm )y ^ j = v anc [ = p$ + j + f } 2 

we obtain the equivalent system (see appendix A): 

(IV 

— = 0 u + f(v,u) 
at 

^ = -bu + g(v,u), (40) 

dt 

where / and g are defined in Appendix A. 

By [11], system (3) has a center manifold of the form u = h(v) and the flow 
in the center manifold ( and therefore in the system ) is given by the equation 

v' = f(v,h{v)) ~ f{v,(j>{v)), 

where h(v) = a$v 2 + a±v 3 + 0(v 4 ), and a,;’s are given in Appendix A. Expanding 
the Taylor series of / we obtain the flow equation 

, b 3 (3m. + b 2 f3 2 m. + b 3 r yCX 2 — b 2 /3p5 + b 3 cx/3m + b 3 p5CX 2 — b 3 f3mcx 2 2 sx 

v = -- v + U(v ) 

b b 

= Hv 2 + 0(v 3 ). (41) 


Therefore the dynamics of solutions near the equilibrium E = (0, 0) is given by 
the quadratic term, whenever this term is not zero. We note that H = 0 if and 
only if 

—/ 3(bm + fim — pS + bam ) 

“ 2 = b{ P S + 7 -pm) ' ( ’ 

Substituting again TZq = 1, expressed as /3m = pS + 7 + hi we obtain H = 0 if 
and only if /3 2 = 3 ( 0 : 2 ). 

If (/I 2 , o 2 ) £ A 3 then H > 0. v' > 0 for v/0. If {^ 2 , 02 ) £ A\ U A 2 then 
H < 0, v' < 0 for v 7 ^ 0. In both cases £ is a saddle point. Moreover, if 
(/? 2 , 02 ) £ A\ U A 2 then H < 0 and v' < 0 for v > 0. Recalling v(t) = I(t) we 
have under this assumption that I'(t) < 0 for / > 0 therefore I(t ) —t 0 + , while 
as v\ = (1,0) is the stable direction of the point E then S(t) —t 0, therefore the 
solutions in the region D approach the equilibrium E as t —> 00. □ 
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Figure 5: Phase plane of the system for 7^o = 1 with (/3 2 , 0 : 2 ) £ A 3 


Example 10. Take the following values for the parameters: p = 0.2, a = 
0.4 ,6 = 0.01,7 = 0.01, a 2 = 10, m = 0.3 ,p = 0.02,6 = 0.2, @ 2 = 0.0498. In 
this case TZq = 1, a® = 1-8876 and g(a 2 ) = 0.4040, hence {@ 2 , 0 : 2 ) £ A 3 . 
By the first case of theorem f the system has a unique endemic equilibrium in 
S 2 = 0.11210, 1 2 = 0.4781. By theorem 9 the equilibrium E is a saddle point, 
see figure 5. 


Example 11. If we take the same values as in the previous example except 
a 2 = 2, then g{a 2 ) = 0.0056 < @ 2 , hence {@ 2 ,oi 2 ) £ A 2 . By theorem f the 
system has no endemic equilibria, and by theorem 9 the point E is a saddle 
point. Moreover, the region D is totally contained in the stable manifold, see 
figure 6 . 


3.2 Stability of endemic equilibria 

The general form of the Jacobian matrix is 

_ @1 (3S 

DF = 


1 T otl 

@1 


-b 


PS 


(1 + al) 2 


— bm + pS 


— pS — 7 — 


P 2 


1 + aJ (1 + al) 2 * ' (1 + a 2 I) 2 , 

Therefore the characteristic equation for an endemic equilibrium is 




@S 


@1 

1 + al 


PS 

(1 + al ) 2 


(1 + al) 2 

— bm + pS 


— pS — 7 — 


@2 


(1 + a 2 iy 


- A 


(43) 


(44) 
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Cursor position: (0.657,0.241) 


Figure 6: Phase plane for TZq = 1 with (fi 2 , ot 2 ) £ A 2 
If we denote by 


C, -.= 
Cs~ 

r>! := 


pi 

l + al 

PS 

(1 + al) 2 

P2 


(1 + a 2 I) 2 ' 

Then the characteristic polynomial is rewritten as 

P( A) = A 2 + WX + U. 

Where: 

W = C 7 + b - C s +pS + 7 + £>/ 

17 = C/7 + Ci Dj — 6Cs + bpS + 67 + bDi + Cibm. 


(45) 

(46) 

(47) 

(48) 

(49) 

(50) 


By proposition Routh Hurwitz criteria for n = 2 if the coefficient W and the 
independent term U are positive then the roots of the characteristic equation 
have negative real part and therefore the endemic equilibrium is locally asymp¬ 
totically stable. Note that whenever the equilibriums are positive, Ci,Cg, Dj 
will be positive as well. Let us analyze the stability according to the value of 
IZo- 

Theorem 12. Whenever the equilibrium, E\ exists it is a saddle and therefore 
unstable. 

Proof. Consider E\ = (Si, Id) and its characteristic polynomial (48). By Routh- 
Hurwitz criterion for quadratic polynomials, its roots have negative real part if 
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and only if U > 0 and W > 0 , where U, W depend on E \. Moreover, when 
U < 0 its roots are both real with different sign and when U > 0 and W < 0 
the roots have positive real part. Computing the value of U and expressing Si 
in terms of I\ we obtain 

h{a\Il + bih + ci) _ hF{I\) _ . > 

(l + a/i)(l + a 2 /i ) 2 (l + cdi)(l + a 2 /i) 2 ’ 1 j 

Where: 

°i = aliPl + bpaS + ba 7 + bm/3 ) = a 2 A > 0, 
b\ = 2 a 2 (/?7 + bpaS + ba^ + bm/3 ) = 2 A > 0, 

Ci = /3/3 2 + bm/3 + bpaS + ba(3 2 + /3q — ba 2 /3 2 + ba 7 = B — a 2 C. (52) 


We are assuming that equilibrium E\ exists and it is positive, and these happens 
(by previous section) when B < 0 and C > 0, so Ci < 0. The sign of U is equal 
to sgn(F(/i)). F{I\) has two roots of the form: 


_ —bi + \Jb\ — 4aiCi 

2ai 

-61 - sfb\ - 4 aiCi 
I * * = ----- 

2di 


(53) 

(54) 


Where b\ — 4aiCi > 0 and therefore I* and I * * are both real values with 
I * * < 0. F(Ii) > 0 for Ii > I* and Ii < I * *, but second condition never 
holds because /1 >0, so F(1 1 ) < 0 for 0 < Ji <7*. 

Computing 7* in terms of A,B,C: 

I* = —— + ~^—J(A 2 -a 2 AB + a 2 2 AC). (55) 

0-2 a 2A v 

Substituting A = E> 2 — 4 AC > 0 


u =-i + ^A]K A, -» AB+ 

= _i + _i_ v , dWTF^ 

—B — \fK 
2A 1- 


(56) 


Therefore U < 0 and the equilibrium E\ is a saddle. 


□ 


Theorem 13. Assume the conditions of theorem 4 for existence and positivity 
of the endemic equilibrium E 2 . If I 2 < I* the equilibrium E 2 is unstable, else if 
I 2 > I* then E 2 is locally asymptotically stable for s > 0 and unstable for s < 0. 
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Where s = m\{—B + \JB 2 — 4AC) A 2Am 2 , 

mi = (r A fi 2 a — (3 2 a 2 A 2Ba 2 )A 2 — a 2 rAC — ABa 2 (ba 2 + 2 r + B 2 a 2 r), 
m 2 = bA 2 — ACa 2 (ba 2 A 2?’) + a 2 rBC , 
r = a(p6 A b A 7 ) A /3. (57) 


Proof. Consider i ? 2 = (S 2 ,I 2 ) be real and positive, and its characteristic poly¬ 
nomial (48). We will have that the equilibrium is unstable when U < 0 and 
locally asymptotically stable when U > 0, W > 0. Following the previous proof 

tj _ h{ail 2 + biI -2 + ci) _ _ I 2 F(I 2 ) _ , , 

(1 A a/ 2 )(l A a 2 I 2 ) 2 (1 + a/ 2 )(l + a 2 I 2 ) 2 ' 1 ’ 

Where ai, 6 i,Ci are the same as in previous theorem. Therefore sgn({7) = 
sgn (F(I 2 )). We have seen that F(I 2 ) has two real roots I* and I * *. Again 
F(I 2 ) > 0 for I 2 > /* and I 2 < I * * (which does not holds because I * * < 0), 
and F(I 2 ) < 0 for 0 < I 2 < I* . So if I 2 < I* the equilibrium E 2 is unstable. 
When I 2 > I* then U > 0 and 


W = (l + alJl + aM^ 2 + <« + £ + cpS) h s 

A a 2 (ba 2 A 2 ap5 + 2ba + 2a^ + 2(3) I 2 2 
+ (a p8 + ba + /3 + a 02 — P 2 a 2 A a 7 + 2 ba 2 ) I 2 A b] 

G{h) 

(1 + a/ 2 )(l + a 2 I 2 ) 2 

By using the division algorithm, 

G(I 2 ) = ( AI 2 A BI 2 A C)P(I 2 ) 

A jic.[(r A (3 2 0i — (3 2 a 2 A 2 Bex 2 ^A 2 — oujrAC — ABoc 2 {bo . 2 A 2r A B 2 ot 2 r^)I 2 
A 2 


(59) 


A bA 2 — ACa 2 (ba 2 A 2 r) A a^rBC], 


— (.AI 2 A BI 2 A C)P(I 2 ) A 


TOi / 2 A m 2 
A 2 ■ 


(60) 


Where P(/ 2 ) is a polynomial in I 2 of degree one. Being I 2 a coordinate of an 
equilibrium then AI 2 A BI 2 A C = 0 and 


G(h) = 


m\I 2 A m 2 
A 2 ' 


Hence sgn(W) = sgn(G(/ 2 )) = sgn( miJ |j~ m2 ) = sgn(mi / 2 Am 2 ). Substituting 
the value of I 2 , 

mi / 2 A m 2 = ^-(-5 A VB 2 - 4 AC) A m 2 . 
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It follows that sgn( 7711 / 2 + 7712 ) = sgn(mi (—B + y/B 2 — 4AC) + 2 Arri 2 ) = sgn(s). 
Therefore E 2 is unstable if s < 0 and locally asymptotically stable if s > 0. 

□ 


4 Hopf bifurcation 

By previous section we know that the system (3) has two positive endemic equi¬ 
libria under the conditions of theorem (4) . Equilibrium E\ is always a saddle, 
so its stability does not change and there is no possibility of a Hopf bifurca¬ 
tion in it. So let us analyse the existence of a Hopf bifurcation of equilibrium 
E 2 = (S‘ 2 , 12 )• Analysing the characteristic equation for E 2 , it has a pair of pure 
imaginary roots if and only if U > 0 and W = 0 . 

Theorem 14. System (3) undergoes a Hopf bifurcation of the endemic equilib¬ 
rium E2 (whenever it exists) if I 2 > I* and s = 0. Moreover, if 02 < 0, there is 
a family of stable periodic orbits of (3) as s decreases from 0; if 0,2 > 0, there 
is a family of unstable periodic orbits of (3) as s increases from 0. 

The characteristical polinomial for E 2 has a pair of pure imaginary roots iff 
U > 0 and W = 0. From the proof of theorem 12 we have that U > 0 if and 
only if one of the conditions (i), (ii) is satisfied . 

Although, sgn(IU) =sgn(s), so W = 0 if and only if s = 0. By first part of 
theorem 3.4.2 of [17] the roots A and A of (48) for E 2 vary smoothly, so we can 
affirm that near s = 0 these roots are still complex conjugate and 


dRe{X(s)) 

ds 


«=o 


d_ 

ds 



l_d ( _1_ \ 

2 ds \2A 3 (1 + aq/i)(l + a 2 /i) 2 / 
4A 3 (1 + a 1 / 1 )(l + a 2 /i) 2 ^ °' 


(61) 


Therefore s = 0 is the Hopf bifurcation point for (3) . 

To analyze the behaviour of the solutions of (3) when s = 0 we make a change 
of coordinates to obtain a new equivalent system to (3) with an equilibrium in 
(0,0) in the x — y plane ( see appendix B ). Under this change the system 
becomes: 


dx _ aux + ai 2 y + c x xy + c 2 y 2 
dt 1 + ol\ y + an /2 

dy _ a 2 ix + a 2 2 j/ + c 3 xy + c 4 xy 2 + c 5 y 2 + c 6 y 3 
dt (1 + a\y + 07 / 2 ) (1 + a 2 y + a 2 I 2 ) 

Where the ay’s and cfs are defined in appendix B. 
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System (62) and (3) are equivalent ( appendix B ), so we can work with 
(62). This system has a pair of pure imaginary eigenvalues if and only if (3) has 
them too. As we said before it happens if and only if any of conditions (i),(ii) 
is satisfied and s = 0. Computing jacobian matrix DF{ 0,0) of (62) 


DF( 0,0) 


a n 

1 + a.\l2 

_021_ 

(1 + CX2I2) (1 + 01/2) 


Ol2 

(1 + aq/ 2 ) 


022 

(1 + CX2I2) (1 + O1/2) 


(63) 


Tr(DF(0,0))=Tr(Df(S 2 ,h)), det(Z?F(0, 0 )) = det(D/(5 2 , h))- 

So condition s = 0 is equivalent to an(l + 02 / 2 ) + 022 = 0 and (i),(ii) are 
equivalent to a 22 a .11 — 012021 > 0 . 

System (62) can be rewritten as 


dx 

dt 

dy 

dt 


anx any 
1 + (X1I2 1 + 011I2 

a 2 ix 

(1 + OL1I2) (1 + 02/2) 


Gi(x,y) 

a 22 y . n , \ 

(l + a 1 / 2 )(l + a 2 / 2 ) Gz( ,2/) ' 


Where G 1; G 2 are defined in appendix B. 

Let A = v/det(.DF(0,0)). We use the change of variable u = 
A(1 + 1 a 1 l/2) + A (i+aij 2 ) > t0 obtain the following equivalent system: 


(64) 

(65) 


X, V = 


Where 

Hi(u,v) 


H 2 {u,v) 


( 0 A \(u\ (Hi(u,v)\ 
V _A oy \y) \H 2 (u,v)J 


( 66 ) 


((—ai 2 Ci + ouc 2 ) u + (—Ac 2 oi/ 2 + Aar 2 ai — Ac 2 ) v) ((A + Aar I 2 ) v — anu) 
ai 2 ((aqA + Aai 2 / 2 ) v + a 12 — aq anu + 0120 : 1 / 2 ) 

(67) 

[(A(l + Oi I 2 )v - an u) (Aiv 2 + A 2 uv + A%v + A^u 2 + A 5 u)] 

( 68 ) 


h(u, v) 


And Ai, ..., A 5 , h(u, v ) are dehned in appendix B. Let 

[(Hi) uv ((Hi) uu + (Hi) vv ) 

~~ {H 2 )uv{{H 2 )uu + {Hljvv) ~ {Hl)uu(H 2 )uu + (Hl)vv{H- 2 )vv]- (69) 


02 — [(//l )uuu 4” (T/'l )uvv 4” ( 7/2 )uuv 4” (dd2)vvv\ 4" ^ 


Where 

and so on (a 2 is explicitly expressed in appendix B) . 

Then by theorem 3.4.2 of [17] if <22 7 ^ 0 then there exist a surface of periodic 
solutions, if d 2 < 0 then these cycles are stable, but if 0,2 > 0 then cycles are 
repelling. 
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5 Discussion 


As we said in the introduction, traditional epidemic models have always stabil¬ 
ity results in terms of TZq, such that we need only reduce TZq < 1 to eradicate 
the disease. However, including the treatment function brings new epidemic 
equilibria that make the dynamics of the model more complicated. Now, let’s 
discuss some control strategies for the infectious disease, analysing the parame¬ 
ters of the treatment function (a 2 ,/3 2 ) and looking for conditions that allow us 
to eliminate the disease. We make this study by cases. 

A first approach is focus on the definition of TZq, we can see that TZq decreases 
when P2 increases, so the first measure suggesting control is a big value for /3 2 . 
But this is not always a good way to proceed. Let us divide our analysis in the 
following cases: 

Case 1: There is no positive endemic equilibrium for TZq < 1. This happens 
when 7^5 < 1 ( by theorem 4 ) or when TZq > 1 and (a 2 ,/3 2 ) £ Ai U A 2 ( 
theorem 4 , number 5 ) . In this case if TZq > 1 there is a unique positive endemic 
equilibrium, therefore there exists a bifurcation at TZq = 1 : from the disease free 
equilibrium, which is globally asymptotically stable for 0 < TZq < 1 (by theorem 
6 ) and a saddle for 7Z 0 = 1 and 02 7 ^ 9(0:2) (theorem 9 ), to the positive 
endemic equilibrium E2 as TZq increase. E2 will be locally asymptotic stable or 
unstable depending on theorem 13 or surrounded by a limit cycle (theorem 14 
) . If conditions for Hopf bifurcation hold then the stability of the limit cycle 
is determined by 02 ; when a 2 < 0 the periodic orbit is stable and therefore E 2 
is unstable, while if <22 > 0 then the periodic orbit is unstable and E 2 is stable 
. In this case the best way to eradicate the disease is finding parameters that 
allow TZq < 1, because then all the infectious states tend to / = 0. 

Case 2: There exist endemic equilibria for TZ 0 < 1. This happens when 
(a 2 ,p 2 ) £ A 3 . The existence of endemic equilibria is determined by the rela¬ 
tionship between TZq and max{Pi, TZq }. Let A(a 2 , P2) = P-a — K(f> G(a 2 , /3 2 ) = 
TZq — Pi, and focus on the implicit curves defined by F = 0 and G = 0. These 
curves divide the domain A 3 in another ones (see figure 7 ): 

A 3 = {(02, &) £ A 3 ,0 < TZq < TZq } 

A 3 = {(02, P2) £ A 3 , TZq > 7^q } 

A3 = { (02 , fa) £ A3,0 < IZo < Pi} 

A 3 = {(^ 2 ,^ 2 ) £ A 3 ,7^o > Pi}- (70) 

If (a 2 , 02 ) £ A 3 n A 3 then there exist two endemic equilibria E\( a saddle ) 
and E 2 ( stable or unstable depending on conditions of theorems 13 and possibly 
with a periodic orbit around (theorems 14 )), but when TZq = 1 one of them 
becomes negative, leaving us with A 2 . In this case TZq < 1 is not a sufficient 
condition to control the disease, because even with TZq < 1 we have endemic 
positive equilibria that could be stable and then the disease will tend to a non 
zero value; also we have the possibility of a periodic solution, or biologically, an 
outbreak that will apparently “ disappear ” but will re-emerge after some time. 
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- Ro=Ro- - Ro=Ro+ . Pl=Ro - Ro=l 

- g - - root=0 _ 


Figure 7: Bifurcation diagram in terms of /3 2 and a 2 . The values of the param¬ 
eters taken are a = 0.4, /? = 0.3, & = 0.2,7 = 0.03,5 = 0.05, p = 0.3, m = 0.3. 
Here TZq < TZq inside the solid curve ( TZq = TZq) and TZq > TZq outside 
it, whenever (^ 2 , 012 ) is in its domain (under the long dashed line “root=0”) 
. 1Z 0 < 1 above the dot-lined line and 1Z 0 > 1 under it; 1Z 0 < Pi above the 
dotted line and TZq > Pi under that one. The areas Hi,H 2 ,A 3 are delimited 
by the dashed line g = 0 and a 2 = 3.8. In this case the endemic equilibria E\ 
and E 2 exists both in the area delimited by the line TZq = 1 and the dotted line 
TZ 0 = Pi , while E 2 exists by itself under the line TZq = 1. 
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The best way in this case is ensuring (a 2 , P 2 ) G (A§ n A|) c because then we 
don’t have endemic equilibria for IZo < 1 and the disease free will be globally 
asymptotically stable. 
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A Computing center manifold 


The Jacobian matrix of system (38) is 

DF(», 0 )=(-; - /Jm - 0 fen + M ). (Tl) 


With eigenvalues Ai = —b and A 2 = 0 and respective eigenvectors v\ = (1,0) 
and V 2 = (— 7 +fe+ bm ; l). Using the eigenvectors to establish a new coordinate 
system we define: 



Under this transformation the system becomes 


du 

dt 


dv 

dt 


= | (7 + flg + bm) Tt I (t) 

dt b 


P [s + mPj 



1 H- ct / 


S + m j I 

1 +al 

dl 


dt 


P 

(^S + rn'j I 


— bS — bml + p61 + (7 + Pz + bm) 


- (pS + j) I - 


1 


P 2 I 

1 + ot 2 I I b 


1 Oc I 


- (pfi + 7 ) 1 - y 


a 2 I' 



(73) 

(74) 


Substituting S = u— C+h2+bm)v ^ j _ y anc j ^ ^ + fj 2 we obtain: 


dv 

dt 

du 

dt 


0 w + f{v, u) 
-bu + g(v,u), 


(75) 


25 



where 


,, s v (—/3 b - /3 ba 2 v) u 
/ ( “-'’) = -(l + a„)(l + a s v)b 

— J 7 -, -77 zr-, -77 ((/3bma 2 + byaa 2 - f3 a 2 p5 + bp5 aa 2 + /3 2 a z m) v 2 

(1 + av) (1 + 0 L 2 v) b 

+ ( bp5 a 2 + f3 bm — (3 bmalpha2 + by a 2 + f3 2 m — /?pS + ba (3 m) v), 

g{u,v) = - —- — --~^[v((mb 2 yaa 2 + 2f3 2 bm 2 a 2 + /3a 2 p 2 S 2 

(1 + av) (1 + a 2 v) b- 

+ (3 3 a 2 m 2 — by p5 a a 2 + by a a 2 f3 m — 2 0 1 a 2 mpS + bpS a a 2 (3 m 

— b 2 ma a 2 (3 — f3 2 bma 2 + b 2 m 2 f3 a 2 + b 2 mp5 a a 2 — bp 2 S 2 a a 2 

— 2(3 bma 2 pS — b 2 m(3 a 2 + bf3 a 2 p5)v 2 + ( b 2 m 2 f3 — 2/3 bmp6 — (3 b 2 m + 2 (3 2 bm 2 

— f3 2 bm + f3p 2 6 2 + 2/3 bma 2 pS — bp5 a/3m + /3 3 m 2 + (3 bp6 — b 2 a (3 m — 2 f3 2 mp5 

+ b 2 m 2 a (3 + ba f3 2 m 2 — uf3 b 2 ma 2 — b/3 2 ua 2 m + b(3 ua 2 pS — y bpS a 2 + y (3 bma 2 
+ b 2 mp6 a 2 — b 2 m 2 /3 a 2 + u/3 b 2 a 2 — bp 2 8 2 a 2 — (3 2 bm 2 a 2 + b 2 my a 2 )v — b 2 m.f3 u 
+ u/3b 2 — bf3 2 um + b(3 up5)]. (76) 

By [11] the system (75) has a center manifold of the form u = h(v). Let 
0 : M —>• M. and define the annihilator: 

N4> = (j>(v))) + b(j>- g(v, </>(v)) 

= -77---77--7 \bpS a v 3 a 2 j3 m + b 2 m 2 f3 v 2 — (3 v 2 b 2 m + b 3 d> + b 3 (h a v 

tr (1 + av)(l + a 2 v ) 

+ b 3 (j>a 2 v + b 2 mya 2 v 2 + <j>(3vb 2 + vb(j)f3p5 + b 2 mp5 a 2 v 2 — tj> f3v 2 b 2 m.a 2 

— 7 bpS a 2 v 2 + b 2 my a v 3 a 2 +7/3 v 2 bma 2 + b 2 m 2 a v 2 f3 — 2 f3 2 v 2 mp8 — (3 2 v 2 bm 2 a 2 
+ /3 v 2 bp8b 2 v 2 a f3m + ba v 2 (3 2 m 2 — bpS a v 2 (3 m + (3 3 v 2 m 2 — 2/3 v 2 bmpS + (3 3 v 3 a 2 m 2 
+ 2 (3 2 v 2 bm 2 + /3 v 2 p 2 5 2 — (3 2 v 2 bm — (f>/3 vb 2 m — vb(j> (3 2 m + (3 v 3 ba 2 p8 — b 2 v 3 a a 2 /3 m 

— bp 2 S 2 a v 3 a 2 — 2 (3 2 v 3 a 2 mp8 + 0/3 v 2 b 2 a 2 — b 2 m 2 (3 v 2 a 2 + b 2 m 2 /3 v 3 a 2 — /3 v 3 b 2 ma 2 

— f3 2 v 3 ba 2 m. + 2 (3 2 v 3 bm 2 a 2 — bp 2 5 2 v 2 a 2 + /3v 3 a 2 p 2 S 2 + b 3 <j>av 2 a 2 — v 2 b(j> f3 2 ma 2 
+ v 2 b(f> (3 pS a 2 + by a v 3 a 2 (3 m + b 2 mp5 a v 3 a 2 — 7 bp5 a v 3 a 2 —2/3 v 3 bma 2 pS 

+ 2/3v 2 bma 2 pS). (77) 

Assume that </> = a^v 2 + a±v 3 + 0(v 4 ), then by substituting <j> and ^ in the 
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annihilator N<j> and expanding its Taylor series we get: 

N(j> = -^-((7/3 bma 2 + b 2 mpS 0:2 — & 2 m 2 /?og + 2/3 iunog pd + 2(3 2 bm 2 + b 2 m 2 (3 
b z 

— ft b 2 m + b 3 ao — (3 2 bm 2 a2 —2/3 bmp5 + b 2 mr)ot 2 — 'ybpd og — b 2 a(3m — (3 2 bm 
+ ba f3 2 m 2 + b 2 m 2 a (3 — 2 (3 2 mpS + (3p 2 5 2 — bpd a (3 m — bp 2 d 2 a 2 + (3 bpS 

+ (3 3 m 2 )v 2 — -^[a(3 3 m 2 — a 0 (3b 2 — b 3 a j — 2 bp5a(3m — (3 2 bm 2 a2 2 — bp 2 d 2 a 2 2 
b z 

+ mb 2 7 a 2 2 — b 2 m 2 (3 ag 2 — b 2 a (3 m — b 2 a 2 (3 m — a (3 2 bm + ba 2 (3 2 m 2 + b 2 m 2 a 2 (3 
+ a (3p 2 5 2 + 3 ao (3 b 2 m, + 3 ao b(3 2 m + 2 a 0 6 2 7 ag +2 ag b 2 pS a .2 — 2 ao (3 b 2 maz 

— 3 ao b(3p5 + 2 ao b 2 a (3 m — 2 a (3 2 mpS + a (3 bpS + b 2 mp5 ctg 2 + 67 ag 2 /3m 

— bjpS a 2 2 + 2/3 bma 2 2 pS — bp5 a 2 (3 m + b 2 m 2 a (3 + 2 ba (3 2 m 2 ]v 3 + 0(i> 4 )). 

(78) 

By choosing the coefficients of v 2 and v 3 in order to have N<j> = 0(v 4 ) we obtain 
that ao and a\ must be the following: 

ao = — — [b 2 m 2 (3 + (3 bpd — b 2 a (3 m + b 2 mpS 0:2 — 7 bpS 02 + 7/3 bma 2 — 2 (3 bmpS 
tr 

— (3 b 2 m + 2 (3 2 bm 2 — (3 2 bm + (3p 2 5 2 + 2/3 bma 2 pS — bpd a(3m + (3 3 m 2 — bp 2 d 2 a 2 

— b 2 m 2 (3a2 — 2 (3 2 mpd + b 2 m"/a 2 + ba(3 2 m 2 + b 2 m 2 a(3 — /3 2 6?n 2 o 2 ], (79) 

ai = p [a / 3 3 m 2 — ao (3 b 2 — 2 bpd a (3 m — (3 2 bm 2 a2 — bp 2 d 2 a 2 2 + mb 2 7 a 2 — & 2 to 2 /3o 2 

— b 2 a (3 m — b 2 a 2 (3m — a (3 2 bm + ba 2 /3 2 m 2 + b 2 m 2 a 2 (3 + a(3p 2 d 2 + 3 ao (3 b 2 m 

+ 3ao b(3 2 m + 2 ao b 2 ja 2 +2 ao b 2 pd 02 — 2 a 0 (3 b 2 ma 2 — 3 ao b(3 pd + 2 ao b 2 a (3 m 

— 2a(3 2 mpd + a/3 bpd + b 2 mpd ag 2 + 67 0g 2 /3m — &7p<5ag 2 + 2 (3bma2 2 pd 

— bpd a 2 (3 m + b 2 m 2 a (3 + 2ba (3 2 m 2 ]. (80) 

Hence h(v) = a 0 v 2 + ait > 3 + 0(v 4 ). 

B Hopf bifurcation 

To analyze the behaviour of the solutions of (3) when s = 0 we make a change 
of coordinates x = S — S 2 , y = I — I 2 , to obtain a new equivalent system to (3) 
with an equilibrium in (0,0) in the x — y plane. Under this change the system 
becomes in: 

dx anx + a\ 2 y + c\xy + C 2 y 2 + C 7 
dt 1 + Oij/ + 01/2 

dy_ _ a 21 X + a 22 y + c 3 xy + c 4 xy 2 + c 5 y 2 + c 6 y 3 + c 8 
dt (1 + oiy + oi/ 2 ) (1 + o 2 j/+ 02 / 2 ) 
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Where: 


an = — b — /?i/ 2 — bot\I 2 (82) 

ai 2 = —2 bma\I 2 + bmai — baiS 2 + 2 p6 a\I 2 + pS — bm — f3\S 2 (83) 

ci = -6ai - /?! (84) 

c 2 = —bma i +p5 oq (85) 

a 2 i = — 1 2 (—/3i — P\Qt 2 I 2 ) (86) 

a 22 = -2p<5«i/ 2 + 2 ( d 1 a 2 £' 2 / 2 - 3 pS oqa 2 / 2 2 - 2 ^ail 2 

— 27 a 2 / 2 — 2pS a 2 I- 2 — 2f3 2 ail 2 — 370 : 1 a 2 / 2 2 — 7 — - (3 2 + (i\S 2 (87) 

C3 = 2 / 3 ia 2 I 2 + (88) 

c 4 = !3ia 2 y 2 (89) 

c 5 = —3p<5 aia 2 / 2 — 37 a! a 2 I 2 — p5 a 1 + (3\a 2 S 2 — 7 ai — j a 2 — pS a 2 

— /3 2 ai (90) 

Cq =—p5 a\a 2 —'y a\a 2 (91) 

C 7 = —(/3iS 2 I 2 - bmail 2 + bS 2 — pS I 2 - pS a\I 2 2 + baiS 2 I 2 + bml 2 

— bm + bma\I 2 2 ) (92) 

c 8 = -I 2 [pS «i / 2 + pS + p6 a 2 I 2 + 7 a 2 I 2 - fd\a 2 S 2 I 2 + 7 ai / 2 + /3 2 otiI 2 

+ 7 + 7 «i a 2 I 2 2 — PiS 2 + f3 2 + pS a\a 2 I 2 2 ]. (93) 


But from the equations for the equilibrium point we can prove that C 7 = c 8 = 0, 
so the system we will work on is 

dx _ ana: + a\ 2 y + C\xy + c 2 y' 2 
dt 1 -\- aiy *b 0 : 1 1 2 

dy _ a 2 ia: + a 22 y + c 3 xy + c^xy 2 + c 5 y 2 + c 6 y 3 , g4 > 

dt (1 + axy + oq/ 2 ) (1 + a 2 y + a 2 I 2 ) 

If we denote system (3) as ( S , /)' = f(S, I ) and system (81) as (x , y)' = F(x, y), 
f = U 1 J 2 ), F = (Fi,F 2 ) then 

F(x, y) = f(x + S 2 ,y + I 2 ), 


and 


0F t 

dx 


(x,y) 


dfi. „ T \ OS, , dfi, „ T \dl, , 

^ (a:+5 2 , y+/ 2 ) ^ (a:, y)+ {x+S 2 , y+/ 2 ) 7 ^ .(*, y) 


dfi 

dS 


(x+S 2 , y+I 2 ) 


m 

dy 


(x, y) 


dfi, „ T ,dS, , dfi, „ r , dl, . 

(x+S 2 , y+I 2 ) — {x,y)+-^Y (x+S 2 , y+.7 2 ) — (a;, y) 


gji 

95 


(a :+S 2 ,y+I 2 ). 


So, the jacobian matrix of (62) DF{ 0,0) in the equilibrium is equal to the 
jacobian matrix of system (3) Df(S\,Ii). We can also compute the partial 
derivatives of system (81) and (62) to prove that they are equal,ie, 


Df(S 2 ,I 2 ) = DF( 0,0). 


(95) 
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Therefore the system (62) and (3) are equivalent and we can work with system 
(62). The jacobian matrix DF( 0,0) of (62) is: 


DF( 0,0) 


a n 

1 + CX 1 I 2 
021 

(1 + a 2 / 2 ) (1 + OC 1 I 2 ) 


Ql2 

(1 + ai/ 2 ) 


022 

(1 + CX 2 I 2 ) (1 + £* 1 ^ 2 ) 


So system (62) can be rewritten as 


(96) 


dx 

dt 

dy 

dt 


a 11 X oi 2 V 

1 + On I 2 1 + U 1 I 2 

CI 21 X 

(1 + OL 1 I 2 ) (1 + CX 2 I 2 ) 


Gi(x,y) 

022 y t , 

' (l + a 1 / 2 )(l + a 2 /2) Gz( ,2/) ' 


Where 


(97) 

(98) 


G\ = - --- . -rrlK 1 + 0:1/2)01 - anai}xy + [c 2 (l + aq/2) - axa^y 2 } 

(1 + aiy + ai/ 2 )(l + ai/ 2 ) 

(99) 

^2 = T—-—- T w, -—- J W, - r w-, -rr{[ c 3(l + «ld2)(l + CX 2 I 2 ) 

(1 + «i y + aii 2 )(l + a 2 2/ + a 2 i2)(l + aiJ 2 )(l + CC2-/2) 

- 02i(«2 + oil + 2aia 2 / 2 )]a;y + [04(1 + ai/ 2 )(l + 0:2/2) — a2iOia 2 \xy 2 

+ [05(1 + aq/2)(l + 0-212) ~ 022(0:2 + Oi + 2aia2di)]y 2 + [cg(l + ai/ 2 )(l + 02I2) 

- a220ia 2 ]y' i }■ (100) 

We need the normal form of the system (62). The eigenvalues of DF( 0,0) 
when S 2 = 0 and (i),(ii) are satished are: 

hi, —hi. 


With complex eigenvector 


-1 


-1 


v= — Ai(l + ail 2 ) + an , v=\ Ai(l + ai/ 2 ) + an 


Ol2 


Ol2 


Using the Jordan Canonical form of matrix DF{ 0, 0) and the procedure in [12] 
(p. 107, 108) we use the change of variable u = x, v = 1 to 

obtain the following equivalent system: 



( 0 A \ M , (Hi {u,v)\ 

V _A \ V J + \ H 2(u, v)J ' 


( 101 ) 
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Where 


((—ai2Ci + anc 2 ) u + (—Ac 2 ai/ 2 + Aai 2 ai — Ac 2 ) v ) ((A + Aai/ 2 ) u — an«) 

ai 2 ((aiA + Aai 2 / 2 ) t> + ai 2 — aifluu + ai 2 Oi/ 2 ) 

(102) 

77-7 [(A(l + ol\I 2 )v - an u) 

h[u,v) 

(103) 

A\ = A 2 (1 + ai/ 2 ) [—ai 2 C60 2 / 2 2 ai — anc 2 ai/ 2 2 a 2 ‘ — diic 2 ttil 2 a 2 

— ai 2 C6<a;i/ 2 + anai 2 aiai 2 “/ 2 + anai 2 aia 2 + ai 2 a 22 0i0 2 

— a iiC 2 o 2 2 / 2 — ai 2 C6<a; 2 / 2 — flnc 2 a 2 — ai 2 Cg] 

A 2 = —A (1 + a\I 2 ) [anai 2 Cia 2 2 ai/ 2 2 + ai 2 2 C40 2 / 2 2 ai — 2 ai 2 anC60 2 / 2 “Oi 

2 2 2 2 2 2 

— 2c 2 qi/ 2 a 2 an + ai 2 aia 2 an J 2 — 2 ai 2 anC60i/ 2 — 2 c 2 Q'iI 2 a 2 an 

+ ai 2 2 C 4 «i/ 2 + anai 2 cia 2 Oi/ 2 + anai 2 Ci0 2 2 / 2 + ai 2 “C 4 a 2 / 2 — 2ai 2 anC60 2 / 2 

2 2 2 2 2 

— 2 an c 2 ai 2 / 2 + ai 2 C 4 — 2 an c 2 a 2 — 2 ai 2 anC 6 + anai 2 cia 2 + ai 2 aia 2 an 

— ai 2 2 aia 2 iai 2 + 2 ai 2 ana 22 Oia 2 ] 

A3 = A (1 + 01J 2 ) ai 2 [—ai 2 0507/2 2 a 2 + ai 2 anaia 2 2 / 2 " + 2 ai 2 a 22 aia 2 / 2 

+ 2 ai 2 anQia 2 I 2 — ai 2 C50i/ 2 + ai 2 a 22 ai + anai 2 «i — ai 2 c.50 2 / 2 + ai 2 a 22 a 2 

— 2 anC 2 ai/ 2 2 o; 2 — anC 2 Oi/ 2 — anc 2 — anC 2 o 2 2 / 2 “ — 2anC 2 o 2 / 2 ] 

A4 = — an [—ai 2 2 C40: 2 / 2 — anai 2 Cia 2 Oi/ 2 + c 2 ai/ 2 a 2 an 2 — ai 2 2 C4a 2 / 2 2 Oi 

2 2 2 2 2 

— ai 2 C4 Q;i/ 2 — ai 2 C4 + ai 2 Oia 2 ia 2 — ai 2 ana 22 aia 2 — anai 2 cio 2 ai/ 2 

+ ai 2 anC 6 + an 2 c 2 a 2 2 / 2 + ai 2 anC60 2 / 2 2 ai + an 2 c 2 a 2 — anai 2 Ci « 2 
+ ai 2 anC 6 Q ! 2 / 2 + ai 2 anC60i / 2 + c 2 ai/ 2 2 a 2 2 an“ — anai 2 Cia 2 2 / 2 ] 

A5 = ai 2 [2 an 2 c 2 a 2 I 2 + an“c 2 ail 2 + ai 2 2 Oia 2 i + an 2 c 2 a 2 2 / 2 2 — ai 2 ana 22 a 2 

— ai 2 ana 22 ai — ai 2 2 C3a 2 / 2 — ai 2 2 C30i/ 2 + anai 2 C5 + ai 2 2 a 2 a 2 i — ai 2 anCi 

+ 2 ai 2 2 aia 2 ia 2 J 2 — ai 2 anCio 2 2 / 2 “ — ai 2 anCiOi/ 2 — 2ai 2 anCio 2 / 2 

+ anai 2 C5a 2 / 2 + anai 2 C50i/ 2 — ai 2 2 C30i/ 2 2 a 2 — ai 2 “C3 + anai 2 C50i/ 2 2 a 2 

2 3 2 2 2 

— 2 ai 2 ana 22 aiQ ! 2 / 2 — 2 ai 2 anCiai / 2 o 2 — ai 2 anCiOi / 2 a 2 + 2 an c 2 Oi / 2 a 2 
+ an 2 c 2 ai/ 2 3 a 2 “ + an 2 c 2 ] 

h(u , v) = A (1 + 0:1/2) oi 2 [(oiA + Aai 2 / 2 ) v + ai 2 — Oian u + ai 2 ai/ 2 ] 

[(a 2 A + a 2 A ai/ 2 ) v + a 12 — a 2 an« + 02/2012] (1 + a 2 / 2 ). 



Hi(;u,v) = - 
H 2 (u,v) = - 

And: 
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Let 


0,1 


10 [(-^L )uuu T (^i)uot T (Hl) UU V T ( Hl)vvv\ 4 " 10^ [(L/l)lil>((L^l)tlU T (-ffl)vtf) 

(H 2 ) UV ((H 2 ) UU + (H 2 ) vv ) - (tfi)„ u (tf 2 ) uu + (Hi) TO (iL 2 )w]- (104) 


Then 

3 ((— CiAvai 2 I 2 + Kvauai 2 — ai 2 ciai/ 2 + auc 2 ail 2 — C\Ava\ - ai 2 Ci + anc 2 ) ai 2 an 2 ai) 

0(2 T q 

8 (ai 2 + aiAii) (1 + ctil 2 y 

(—3 anc 2 — 3 CLnc 2 otiI 2 + 2 a\ 2 a\\cx\ + ai 2 Ci + ai 2 Ciai/ 2 ) a iA 2 
8 (1 + ol\I 2 ) a\ 2 3 

- -1----o [2 anA^ai A + 6 anA^a iA a 2 / 2 + 2 anA5ai 2 A J 2 

8 A (1 + ail 2 ) 4 oi 2 4 (1 + a 2 I 2 ) 3 

+ 4 an A 5 CI 1 2 A / 2 2 a 2 + 2 an A^a 2 A — an” ^ 30:1 — 2 an" A. 3 aia 2 / 2 — an 2 A^a 2 — an A 2 ai 2 
— anA 2 ai 2 a 2 I 2 — ai\A 2 ai 2 ot\I 2 — anA 2 ai 2 aiI 2 2 a 2 + A^Acin + A±Aai 2 a 2 I 2 
+ 2 A4A ai 2 ai/ 2 + 2 A4A ai 2 ai/ 2 2 a 2 + A4A ai 2 ai 2 / 2 " + A4A ai 2 ai 2 / 2 3 a 2 ] 

3 (—Aiai 2 — Aiai 2 a 2 I 2 + A 3 CI 1 A + 2 A^ol\Acx 2 I 2 + A^a 2 A) 

8 (1 + ai J 2 ) 2 ai 2 4 (1 + a 2 / 2 ) 3 

1 r ^ A (—2anc 2 — 2 anC 2 ai/ 2 + ai 2 ano:i + ai 2 Ci + fli 2 Ciai/ 2 ) 

16A ai 2 4 (1 + ai / 2 ) 2 

2 (A 5 A + A^Aa\I 2 — an^) (—anA .5 + A 3 A + AzAaiI 2 ) 

A 2 (1 + ai / 2 ) 6 ai 2 6 (1 + a 2 / 2 ) 4 

^ (—ai 2 ci + anc 2 ) an 2 A .5 (—c 2 «i / 2 + ai 2 ai — c 2 ) A 2 A 3 . . 

ai 2 5 (1 + o;i / 2 ) 4 A (1 + a 2 / 2 ) 2 ai 2 5 (1 + ai / 2 ) 2 (1 + a 2 J 2 ) 2 


31 



